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ABSTRACT 


There are a number of military applications in which the geographic location of a 
signal of interest is of prime importance to the ability of a unit to fulfill its mission. The 
accuracy of the geographic fix provided to the warfigjiter can directly affect the success or 
failure of a mission. One method to improve the accuracy of existing systems is to use the 
wei^ted average of a number of intercepts. Each intercept is manifested as an error ellipse 
comprised of a latitude, longitude, semi-axes, heading and a related Chi-squared distributed 
probability. Individual error ellipses can be viewed as a quadratic surface perpendicular to 
the x,y plane of a bivariate normal distribution, the z-axis intersection of which 
corresponds to a Chi-squared value. By transforming the individual error ellipses to their 
related location covariance matrices, Gaussian statistics may be used to produce a single 
location ellipse that combines information from several two-dimensional target location 
ellipses. By providing a means to fuse data from a given source the warfi^ter or analyst will 
be able to more accurately assess a threat and respond. 
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I. INTRODUCTION 


A. JUSTIFICATION 

There are a number of military applications in which the geographic location of a 
signal of interest is of prime importance to a unit*s ability to fulfill its mission. The accuracy 
of the geographic fix provided to the warfi^ter can directly affect the success or failure of a 
mission. It is for this reason we must continually strive to improve our ability to generate 
the most accurate and condensed geographic solution from all available data. 

This thesis describes an algorithm to produce a single location ellipse that is the 
combination of several two-dimensional target location ellipses. The combination algorithm 
presented in this document is limited to a stationary single source for data set production. 
The purpose of this study is to assist the Naval, Command, Control and Ocean Surveillance 
Center Research, Development, Test and Evaluation Division’s (NRaD), Integrated Satellite 
& Link Communications Division (Code 841) in the documentation of the combination 
algorithm. This effort is in support of the Classic Crystal project developed by NRaD. 

B. SCOPE 

The scope of this thesis involves the development of die combination algorithm. A 
mathematically rigorous presentation of assumptions made in the algorithm and of the 
supporting theory will be presented. If an appropriate method is available, verification of 
developed mathematical equations is supplied. It is not within the scope of this project to 
consider the effect of system dependent biases (e.g. atmospheric delays, refraction) on error 
ellipses produced by different systems. For this reason the data for this project is limited to 
error ellipses produced by a single source. 


I 
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C. CHAPTER OUTUNE 

Chapter II presents the underlying statistical theory upon which the Combination 
Algorithm is based. The chapter culminates with the formula for the optimal estimate of 
the target location as presented by J.A. Roecker, who based his work on that of Nelson 
Blachman. The concepts of a bivariate normal distribution and application of Bayes’ 
Theorem are covered as they directly support the formulation of the optimal estimate. 
Chapter III concerns itself with a number of topics all related to applying the optimal 
estimate of the target location developed in Chapter II to a real world scenario. 
Development of the covariance matrices, the effects of adding an observation to an existing 
solution, ellipse rotation and heading issues are covered in rigorous detail. In Chapter IV 
the theory and concepts of the previous two chapters are put to use in the presentation of 
the Combination Algorithm. The process is detailed in a step by step manner that can 
easily be applied to a number of platforms. Finally, Chapter V offers conclusions and 
recommendations for further study. 
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II. STATISTICAL THEORY 


A. RELATIONSHIP OF THE ERROR ELLIPSE TO A BIVARIATE 

NORMAL DISTRIBUTION 


Each of the target location error ellipses to be combbed to produce a sb^e error 
ellipse will be defined by the followbg parameters: center position, headbg and magnitude 
of semi-axes and probability level. These ellipses represent an area that has the specified 
probability that the target lie withb it’s bounds. We will begb by definbg a jobt normal 


distribution of two variables 






as 


<^(x) = k* expj- ^(x - h)’^b(x - n) 


( 2 . 1 ) 


where k is determbed by the normalization of the total probability to one, H and x are 2- 


component Column vectors and B a (2x2)-matrix. The vector |X - 
of the center of symmetry of the distribution 




gives the location 


J |(x-p)<^x>&,t&2 =0 


( 2 . 2 ) 


The expected or mean value of a random vector x is 


E(x) = I I x^x)cfeit36r2 = 


\l^) 


(2.3) 


The covariance matrix of the distribution is defined as 

C(x) = E(x - pX* - = [e(^, - m)(^; - Mj)] (2-4) 
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(Anderson,1984). The i th diagonal element of this matrix, E^x,. —//, ) , is the variance of 
and the i, j th off-diagonal element, e(x,. - )^x^ - , is the covariance of x, and Xy, 

Note that since e(x,.-/ 4 )^Xy-//yj = E(xy-//yj(x;-/ 4 ) the covariance matrix is 
symmetric. With this said the covariance matrix can be written as 


C = B ‘ 


1(^2-aX^1-a) (^2-aX J 



/ 2 



( 2 




<^12! 



<^n 


\^ 2 t 

<r;j 


^*^12 

gIj 


(2.5) 


( 2 . 6 ) 


where we have used the symmetry of the covariance matrix (i.e. t^^trix 

inversion we obtain 



(7^ (J 2 


i_ [ <^2 -^u 


(2.7) 


To facilitate the next step in our development of the ellipse of covariance we will introduce 
two reduced variables 


/-u 

(^4 


( 2 . 8 ) 


with the property var(M,) = var(« 2 ) = 1» and the correlation coefficient 


p = ^ = coy{u„u,). 


(^ 1^2 


(2.9) 


Equation 2A now takes the reduced form 


,"2) = ^ * expf - 


( 2 . 10 ) 


4 



with 


B = 


l-p 


J ^ 


( 2 . 11 ) 


Lines of constant probability density are determined by requiring the exponent in Equation 
2.10 to be constant, as shown in the following equation: 


+ «2 - 2m,W 2/7) = c 


( 2 . 12 ) 


If we choose the constant, c=l, then in terms of the original variables, Equation 2.12 
becomes 




cr. 


(2.13) 


(Brandt, 1983). This can also be written as a sum of squares of two stochastically 
independent variables 


1 



( A 

iW ^ 

l-p' 

ts) 

1 

t; 

1_ 



+ 


(2.14) 


which are distributed as with two degrees of freedom. This is the equation of an ellipse 
with the center located at p, pz- The semi-axes of the ellipse have an angle 9 with respect 
to the Xj, Xj axes. This angle, and the magnitude of ihe semi-major axis a and the semi¬ 
minor axis bean be derived from Equation 2.13 using the properties of conic sections: 


tan2^ = 


2pa^ 



(2.15) 


_ 0 - 1^2 _ 

crl cos^ 9-2f\(ji sin 0 cos ^ -I- of sin^ 9 


(2.16) 
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(2.17) 


_ 

al sin^ 9 - sin Ozo^d + of cos^ 9 

The ellipse we have just described is known as the ellipse of covariance of the 
bivariate normal distribution or, in the context of this thesis, the error ellipse. The ellipse 
will always fall within the rectangle defined by points (///, and the standard deviations 
(®i. ^ 2 )- It can been seen from Figure 2.1 that the ellipse will touch the rectangle at four 
points and in the extreme case of p = +1 the ellipse will degenerate into a straight line along 
one of the diagonals of the rectangle. To help visualize the concept, the bivariate normal 
distribution density function can be thought of as a surface above the plane. The contours 
of equal density are contours of equal altitude (equal probability) on a topographical map; 
they indicate the shape of the hill (or probability surface) (Anderson, 1984). 



Figure 2.1 COVARIANCE ELLIPSES 


The form that the ellipse of covariance was presented in Equation 2.14 did well to 
illustrate the properties of an error ellipse. In order to implement the combination 
algorithm we must be able to derive the covariance matrix from the initial parameters. This 
is most easily accomplished by rewriting Equation 2.14 in matrix form as 

X^={^- C“’(x - x) (2.18) 
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where x is a Gaussian random vector, x is the true mean of the vector x and C is the 
covariance matrix of x (Roecker, 91). The chi-squared value with 2 degrees of freedom 
corresponds to the probability that the observation falls within the ellipse described by this 
equation. In Equations 2.15 thru 2.17 the angle between the ellipse axes and Xi, X 2 axes, 
and the semi-axes were derived from Equation 2.14. The same information can be 
extracted from Equation 2.18 with the help of the eigenvalues and eigenvectors of C 
(Roecker, 1991). The eigenvalues are found from 

lC-el| = 0 (2.19) 

where I is the identity matrix. Substituting for C from Equation 2.6 and taking the 
determinant yields 


e^-s(c7^+ -4 = 0. 

(2.20) 

Therefore, the eigenvalues are 


gI + al ± 

2 

(2.21) 

and the semi-axes are the square roots of the eigenvalues multiplied by the 
freedom x^ -value corresponding to the specified probability p 

two-degree of 


(2.22) 


(2.23) 


Values of x^ vs. p are given in Table 2.1 
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Table 2.1 Values of Chi-square vs. Probability 



P 

2.30 

0.683 

4.61 

0.90 

6.17 

0.954 

9.21 

0.99 

11.9 

0.9923 

18.4 

0.9999 


B. BAYES* THEORM 

This project is intended to combine target error ellipses (discussed in the previous 
section) generated by several independent observations, each of which can be considered 
Gaussian. The result of which will be an error ellipse that represents the smallest area 
containing the target location with the specified probability p. The nucleus of the 
Combination Algorithm is based on the work of Nelson M. Blachman of GTE 
Government Systems Corporation. Since Blachman’s application of Bayes’ Theorem is 
fundamental to the solution, excerpts of his manuscript (Blachman, 1989) are presented 
below: 


We suppose that, on the basis of I independent sets of observations Oj, O 2 ,..., Oj of 
the same target, a location ellipse has been found. Each such ellipse is a contour of the 

Gaussian conditional probability density function (pdf j of the target 

location based on one set O, of observations. Each of the I given ellipses is based on 
the Bayesian relationship 
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(2.24) 



Pi(x,y)i 

Jo, 

<x,y) 

A 

0.) 



expressing the probability density function (pdf) of the target conditioned on the 

observations 0^ as the ratio of the product of the prior (i.e. before O^) pdf of 
[x,y) times the pdf of O, conditioned on the location [x,y)(this product being the 
Joint pdf of ofid Oj to a normalizing constant, which is the unconditioned pdf 

of Of and is the integral of the foregoing product over all \ and y. Because different 
observers may utilize different prior pdfs, can depend on i. As long as 

p.{x,y^ is constant over the entire area where ^0, it cancels out here 

because is proportional to it, and so its value is not important. 


For each I the pdf of O, conditioned on target location 



is therefore 



p{o)i{x,y\o) 

Pi{^,y) 


(2.25) 


On account of the statistical independence of the different sets of observations for a 
given target location, the product of p[o\ over all i is the conditional pdf 

p{p^, 02 ,..., 0 \x,)^ of the entire set of observations for a given target location 
Multiplying it by the prior pdf p{x,y^, we get the Joint pdf 
p{p^,02,...yOj ,X,J^ of the observations and the target location. Dividing by 
p{p^, 02 , ..,0l), which is the integral of the latter over all x and y, we get, by 
Bayes’ theorem, thepdfof{x,y) conditioned on aH observations. 


9 



(2.26) 


„LJo o] _ £(£^,n£(2LkhM 


The factors M) and which depend only on the 

observations, serve merely to normalize this posterior pdf; they do not effect its shape. 
As long as the prior pdf p{x^y^ is constant over the whole area where the target may 

lie, it cancels out as before, and it too does not ajfect the shape of 

which is thus determined as Equation 2.26 indicates, by the product of the pdfs 

p{^yy\Oi ) based on the individual sets of observations. 


Blachman has set the theoretical foundation from which a practical combination 
algorithm can be developed. In order to apply Equation 2.26 to a spherical earth we must 
have the ability to manipulate the observations by a constant scaling factor. Roecker while 
expanding on Blachman’s work (Roecker, 1991) noted that, 


as long as the prior pdfs p{X) and p^ (X) are constant over the areas of interest, 
they will cancel out with and p{p^ to affect only 

as a scaling factor. This leaves the conditional pdf completely 

dependent on the product of ^ for its shape, and can be expressed as 


p(x\o„...A)=Ktlp{x(o,) 

i=l 


(2.27) 
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where K is the scaling factor and A' is a vector containing the target location. This 
observation will later allow the observation covariance matrices to be scaled, compensating 
for a spherical earth. 


C. METHODOLOGY 

In order to apply Equation 2.27 a method of extracting an optimal estimate of the 
target location must be developed. Since />(x|o;) is Gaussian, the conditional pdf of the 
target location can be written as 


/?(x|o,,..., 0,,) = a:, n expj- ^(x - o, C:‘ (x - o, )| (2.28) 

/?(xlo„..., 0 ^^) = ATj expj- ®r‘(x - o,) I (2-29) 


where Ki is the scaling factor, x the vector containing the target location, o,. the i th 

observation of x and C the covariance matrix of x (Roecker, 1991). The optimal estimate 
for any symmetric cost function is the same for Gaussian distributions (Van Trees, 1968) 
and results by minimizing the pdf detailed in Equation 2.29. The minimum of Equation 
2.29 is found by setting the first derivative to zero as presented below 


dx 




1 ^0 

2 dx ) 


(2.30) 


where E^ is defined as 


1=1 


(2.31) 


and where 


11 



dx 


dxy 

IEm. 


(2.32) 


Note that M is the number of observations and n is the number of components of the 
vectors x and O (n=2 for this problem). To begin with we will examine the first observation 

0l 




X - 
: 


X,-o„ 

* 

1 


J^n-Oxn_ 


n n 

X, -o„ 

Z5,,(x,-0,,) - 

• 

_ i=\ »=1 

Xn-Ox„_ 


j^\ 1=1 


then by the product rule we have 


1=1 J =1 


since all summing is to n we can make a change of indices for readability 




'-t 1=1 


j=l 


Because the covariance matrix is symmetric we have 


(2.33) 


(2.34) 


(2.35) 


(2.36) 


(2.37) 


(2.38) 
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(2.39) 


dE, 

dx, 




n 

/=1 




) 


= 2^B^{x,-o^,) 

(2.40) 

In matrix notation this becomes 



(2.41) 

. . 

Setting the derivative ^ to zero we have 


2 C”'(x-o,) = 0 

(2.42) 

Ci‘x-C>i =0 

(2.43) 

X= 0^ 

(2.44) 


So the initial estimate for x is the first observation o,, which was to be expected 
Now that the process for a single observation has been established, we can generalize to 
multiple observations. To begin with we expand the reduction variable E to include M 
observations 

Em ( 2 . 45 ) 

m=ivy=l i=l / 

M 

p.46) 


and once again set the derivative to zero 
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/9/r ^ 

/M=l 

M 

2[C;'(l-o.)] = 0 (2.48) 

W=1 

f M \ M 

Z«;'*-Z(«>.)=0 (2-4!') 

^m=l m-\ 

M 

= 0 (2.50) 

m-l 

Which, with a little algebraic manipulation yields 

M f M 

x = Cf](C>„) where C= X®;’ (2.51) 

m=l ^««=1 2 


As depicted in Figure 2.2, x is the vector representing the location of the weighted 
average of M number of o observations and C is the resultant covariance matrix. This is 
die optimal estimate of the target location we desire and the premise upon which the 
Combination Algorithm is based. 



Figure 2.2 OPTIMAL ESTIMATE OF TARGET LOCATION 
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III. CONCEPT DEVELOPMENT 


A. COVARIANCE MATRICES 

1, Scaled Covariance Matrix 

Chapter II laid the theoretical groundwork for the Combination Algorithm. We 
have seen that the target error ellipses we are attempting to combme can be modeled by an 
ellipse of covariance. The covariance matrix being the structure defining the target error 
elUpse. In the real world problem, the error ellipse will have been generated upon the 
spherical earth. In order to reduce die overall problem to a two dimensional one, each 
ellipse will be projected into a plane. Therefore, the covariance matrices will need to be 
modified, compensating for this projection from the spherical earth to the plane. This 
modification is to compensate for the inequality between latitude error and longitude error 
as the latitude increases from zero. The geometry of the projection from a spherical earth 
to a tangent plane is depicted in Figure 3.1. We need to find the distance apia„g (length of 

projected ellipse axes) when given the arc-length along the surface of the sphere sphere 

ellipse axes), where te = radius of earth and 

^sphere = 

(3.2) 


Observing the geometry 


tan« = ^ 

re 

(3.3) 

a plane 

(3.4) 


which leads to the transformation 
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Figure 3.1 PROJECTION GEOMETRY 


The congruence transformation 

C = S^CS (3.6) 

produces the scaled covariance matrix C that has maintained the symmetry of the original 
covariance matrix.. The scaling matrix S is obtained by observing from Figure 3.2 below 
that the arc-length measured along the surface of the earth at constant latitude ^ 

corresponding to a change in longitude AZ is given by cos(^)aA ; and the arc-length 
measured along the surface of the earth at constant longitude A corresponding to a change 
in latitude is given by r^A^. To put latitude and longitude on equal footing (in terms 
of arc-length) changes in latitude are multiplied by and changes in longitude are 

multiplied by cos(^) . Therefore, 



0 ^ 
cos^j 


(3.7) 
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A. 1 


\ ^ 


- 




j \ 



Te 


Figure 3.2 GEOMETRY OF SCALING MATRIX 


Expanding Equation 3.6 we find that 


'a 0 


0 ^ 

lo r^cosijfAcT;,^ 

lo 

cos^J 


(r. 

0 ^ 

^r^(T;^^cos<j> r^<T^cos<f>j 

u 

cos (f>j 


(3.8) 


(3.9) 




(3.10) 


This is the scaled covariance matrix, where (j) = latitude, X = longitude and is the radius of 
the earth. In the section to follow it will be necessary to solve for the eigenvectors of this 
matrix. To aid in that process, the following notation will be used to refer to the scaled 
covariance matrix in Equation 3.10. 


y\ 

P) 


(3.11) 
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where use has been made of the fact that cr^ = • 


2. 2x2 to 3x3 Covariance Matrix Conversion 

The objective is to convert the 2x2 covariance matrix given in terms of latitude and 
longitude (Equation 3.11) to a 3x3 covariance matrix in terms of (x,y,z). The transformation 
is as follows 




(3.12) 


where the 3x3 covariance matrices are 





\ 












^zzJ 


r 


\ 



^ijiX 








<=^hX 

^hhJ 


(3.13) 


(3.14) 


and the 3x3 transformation matrix T is given by 


T = 


^ dx 

dx 

dx'" 

at, 

dX 

dh 

dy 

dy 

dy 

at, 

dX 

dh 

dz 

dz 

dz 



Jhj 


(3.15) 


where for a spherical earth 


x = (r,+^)cos^ cos A 
^ +a)cos^ sin/l 


(3.16) 

(3.17) 
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z = +h)sin^ 


(3.18) 


Te being the radius of the earth and the partial derivatives are 


d(f> 


-(r^ + h) sin(^^ )cos(A) 


dx 

lx 


= -(r^ + h) cos(^ )sin(2) 


dx 

'dh 


= cos(^)cos(2) 


d(i) 


- sin(^ )sin(A) 


dx 


= (r^ + A) cos(^ )cos(A) 


dy 

— = cos(^)sin(l) 


dz 

— = (r, +;i)cos(<^) 



dz 

- = sin(«. 


(3.19) 

(3.20) 

(3.21) 

(3.22) 

(3.23) 

(3.24) 

(3.25) 

(3.26) 

(3.27) 


To define we start by finding the unsealed covariance matrix from 

Equation 3.6 


«(„, =(s-Tes- 


(3.28) 
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where 




(I 

re 

0 


0 

1 

r, cos^j 


(3.29) 


and the scaled covariance matrix C is given in Equation 3.11. Using the unsealed 
covariance matrix in terms of latitude and longitude and assuming the heig)it h to be 

zero with zero uncertainty, we can form the matrix 




^ 0 




'JU 


0 oj 


(3.30) 


3. Selection of Eigenvalues and Eigenvectors 

The key factors in the transformation of an observation seen in terms of semi-axes 
a, b and heading 0, to it’s description by a covariance matrix, are it’s eigenvalues and 
eigenvectors. This section will accomplish four tasks: derive the eigenvalues from the 
characteristic polynomial of Equation 3.11, derive various identities from the characteristic 
polynomial that will be useful in later sections, normalize and solve for the eigenvectors, and 
determine which of the solutions will be numerically stable. 

a. Eigenvalues 

In order to produce the eigenvalues, which will be denoted by e, the 
nonlinear equation 


(C-ffI)x = 0 (3.31) 

needs to be solved. To begin the process the scaled covariance matrix C is shifted by dl, 
where I signifies the identity matrix 
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By definition (Strang, 1988) the number £ is an eigenvalue of C if and only if 
det(C- £l) = 0. This determinant leads to the characteristic polynomial. It’s roots, i.e. where 
the determinant is zero, are the eigenvalues. 

det(C - £ l) = (a - £)(y0 - e) - r' 

= £^-(a + y0)£ + ay0-r'=0 (3.34) 


The real solutions of the characteristic polynomial (the eigenvalues) are found by the 
quadratic formula as show below. 



(a + ^)±Ja^ +2aP + |3^ -4ap + 4Y^ 
2 


(3.36) 



b. Identities 

With the characteristic polynomial at hand, this is a good time to observe 
two identities that will become useful in future sections. Both of the following identities are 
obtained from factorization of the general quadratic equation, which is 
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(3.39) 


- [s^ + £L)f + e^s_ = 0. 


(3.40) 


Comparing Equation 3.40 with the characteristic polynomial Equation 3.34 we find two 
identities that will later aid in establishing relationships between the scaled covariance matrix 
Equation 3.11 and the primary data parameters a, b, and 0, 


s. + s_ = a + ^ 


(3.41) 


e^s_ =aj3-/^ 


(3.42) 


Also worth noting, but of less significant impact, are the identities 




(3.43) 


6^ -a and s_\ .= P- 
+ j-=o ir=o 


(3.44) 


X,- = 


>'x ^ 

•*•11 


mentioned yields 


c. Eigenvectors 

This section will present solutions for the eigenvectors x, where 
, Writing Equation 3.31 in terms of Equation 3.32 and the notation just 


a-s Y 

y r P-e) 






VoJ 


(3.45) 


Equation 3.45 represents a system of two equations with two unknowns. By adding the 
formula 


x^,+xl=l 


(3.46) 
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to the system of equations, we can simultaneously nonnalizing the eigenvector to unit 
length and supply an additional formula that will be used to produce four valid normalized 
solutions. Each pair of equations will generate solutions for the positive and negative roots 
of Equation 3.38. Each solution will be checked for computational stability. The most 
stable vector from each root will eventually be used to populate the eigenvector matrix. 

(1) Solution for . The first solution presented is for x+^, 
which corresponds to the solution of the first equation in Equation 3.45. When 

(a - )x, + ^ 



(3.48) 


The second equation needed to solve this system is supplied by normalizing the 
eigenvectors to unit length using Equation 3.46. By substitution we have 



(3.49) 


Solving for we get 



Substituting Equation 3.51 back into Equation 3.46 yields 



(3.50) 


(3.51) 


(3.52) 
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Therefore, 


Check: 


xi’> = 


yjia-s^y +; 


^ -r ^ 

va - sj 


Cx(’> = 




1 




1 





1 


^ a xV -y ^ 

Y pAa-sJ 


^- ay + ay - ye+'^ 


+ocp- psj 


^ -ys^ ^ 

U+£- -psj 


r \ 




+ y^ 


U_ -p. 


Now Equation 3.41 states 


s_-P = a-e^ , 


therefore, 

and Equation 3.31 holds true. 

In checking for numerical stability note that when y->0 we obtain 
of x+^, therefore x+^ is not numerically stable as y-^0. 


(3.53) 


(3.54) 


(3.55) 


(3.56) 


(3.57) 


for the components 
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(2) Solution for xf ^. The second solution evaluated is xV > 
which corresponds to the solution of the second equation in Equation3.45 when 


=0 

(3.58) 

r 

(3.59) 

The second equation needed to solve this system, as in the previous solution, is 
normalizing the eigenvectors to unit length using Equation 3.46. By substitution 

supplied by 

we see 

+il= = i 
r 

(3.60) 

solving for x^ we get 



(3.61) 


(3.62) 

Substituting Equation 3.62 back into Equation 3.46 yields 


S.-P 

(3.63) 

Therefore 



(3.64) 
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Check: 


Cx 


( 2 ) _ . 


1 fa 




K r j 


as^ -aP + y 


2^ 


Ys.-^r+PrJ 




+r 


ae^ -s^s 


\ rs. J 




+r 


^a-s^ 


\ Y J 


Using the Equation 3.41 


(3.65) 


(3.66) 


(3.67) 


(3.68) 


a-6_ = s^-P 


yields 



X 


( 2 ) 


and Equation 3.31 holds true. 


In checking for numerical stability we note that when y—>0 


X 


( 2 ) 


1 

r-(P-a)1 


'^1 

- 




a-p 

1 0 > 


.0; 


(3.69) 


Hence, "ifY is a numerically stable solution. Therefore we will use Equation 3.64 for the 
solution from the eigenvalue . The next step is to determine which of the two equations 
from the eigenvalue S _. will be used. 
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(3) Solution for xi*’. The next solution presented is for 
which corresponds to the solution of the first equation of Equation 3.45 with £ = 

{a-s_)x^ +7X2 = 0 

(3.70) 

r 

Xj - - " ^2 

a-e_ 

(3.71) 

The second equation needed to solve this system of equations is supplied by normalizing 
the eigenvectors to unit length using Equation 3.46. By substitution we see 

+ \ X2-I 

J 

(3.72) 

Solving for X 2 we get 

"^~(a-e_y+r^ 

(3.73) 

(a-s_) 

^(a-s^y +r^ 

(3.74) 

Substituting Equation 3.74 back into Equation 3.46 yields 


-r 

yl(a-e_y+r^ 

(3.75) 

Therefore, 


1 ( -r) 

(1) _ _ 

■ ia-s.y 

(3.76) 
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Check: 


Cx 


( 1 ) _ 


Now Equation 3.41 states 


^j(a-s_} 

+y 

1 

/ 

^(a-s_y 

2 ' 
+y 

1 


^j(a-£_y 

+y' 

s_ 


)' 

' +y^ 




;'V -y 

Y B)\a-s_J 


-ay + ay -ys_ 

.,2 1 Pr, 




-y 


Therefore, 


Cx<‘^ = 


(377) 


(3.78) 


(3.79) 


(3.80) 


showing that Equation 3.31 holds true. 

In checking for numerical stability we note that when y— >0 




1 


0 ^ fo^ 


a - (3 va - pj 




(3.81) 


producing a stable solution. 
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(2) 


(4) Solution for . The last solution presented is for xl 
which corresponds to the solution of the second equation in Equation 3.45 when s — £_ 


=0 


(3.82) 


X, = 


P-S- 

y 


(3.83) 


The second equation needed to solve this system, as in the previous solutions, is supplied by 
normalizing the eigenvectors to unit length using Equation 3.46. By substitution we see 




y 


^2 =1 


(3.84) 


Solving for Xj we get 


^2 = 


[p-s_) 


(3.85) 




+y 


(3.86) 


Substituting Equation 3.86 back into Equation 3.46 yields 


X, = 


s.-P 


yl{p-e-f+y' 


(3.87) 


Therefore, 


X 


( 2 ) _ 


1 

^(p-s_y ^y"- 



(3.88) 
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Check: 



1 


fa 



-p\ 



2 

1 


P) 

V ; 

y J 


1 


^a£_ 

— ocfi -f 




2 7 

I +y^ 


-Py + 

Pr) 


1 


'as_ 

-8 



a/(a- 

^-y 

' 2 
+r 


rs- 

J 



e_ 

j 

r 

a - 





^y 

' +r^ 

^ r 

J 




(3.89) 


(3.90) 


(3.91) 


(3.92) 


Using Equation 3.41 


a-8^ = s_- P 


in Equation 3.92 gives 


showing that Equation 3.31 holds true. 

Checking for numerical stability we note that when y->0 we obtain for the components 
of , therefore, x ^2) jg ^ot numerically stable as y—>0. 

The eigenvalue£■_ like the eigenvalues^, produced one computationally stable 
eigenvector. Note that s_ > 0 ands^. > 0 so both values are positive. The equations to be 
used in solving for the eigenvector matrix G will be: Equation 3.64 which corresponded to 
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the solution of the second equation in Equation3.45 when £ = , and Equation 3.76 which 

corresponds to the solution of the first equation of Equation 3.45 with e-s_. 


d. Eigenvector Matrix 

This section will develop the solution for eigenvector matrix G and verify it 
by producing the diagonalized eigenvalue matrix A. The eigenvector matrix G is the matrix 

whose columns are the eigenvectors and (Recall these are the computationally 
stable eigenvector solutions found in sub-section c. of this chapter). In order to simplify 
the construction of G, some algebraic manipulation of the vectors is required. Consider 
Equation 3.64. Since - P = a- s_ (Identity Equation 3.41) we have 


1 1 


yl(a-s_)(s^ -p) + r^ 

\ r ) 


(3.93) 


Now consider Equation 3.76. Again, since a — S_ - s^- P we can say that 


r(V _ 


1 

f -r ) 

^{a-s_)(s^-p} + r^ 

la - s_J 


(3.94) 


By combining Equation 3.93 and Equation 3.94 the eigenvector matrix G is 


G = 


yl{a-s_){e^ -p) + f 


s^-p -r ^ 
K r ci-e_y 


(3.95) 


The determinant (difference of the products of the numbers in the two diagonals) of G is 

1 


det(G) = 


^{a-e,)[s^-p) + y^ 
which is equal to one, and the inverse of G is 


\8^-p)[a-e_) + Y^}, ( 3 - 96 ) 
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(3.97) 


G ’ = 


yl(a-s_){e^ -fi) + r 


fa-e_ Y ^ 

rL -y -pj 


The diagonalized eigenvalue matrix A(the diagonal matrix with the eigenvalues s_ and as 
the diagonal components) is 


(3.98) 



\ 

\ 

(a-s_ 

r ^ 


y a-sj\-y p) 

1 -r 

-Pj 


^a-s_ r Yae^-afi + r^ -ay + ay-ys. 



1 

f 

(a-e\s^-p) + y^ ^ 

> f 

{a-s_] 

> f 


a-s_ y 


'^^ae^-E^e_ -ys_ 


a-s_ y I 

Y sX^^-PV 

eXa-e.y +r^e+ -r<e’-(«-'£'-) + ?«'-(^+ “^) 


(3.99) 


(3.100) 


( 3 . 101 ) 


(3102) 


and by using a- e_ = e^- p (Identity Equation 3.41) this becomes 


A = 


X 0^ 

^0 S.J 


(3103) 


With this check it is reasonable to assume that Equation 3.95 is in fact the eigenvector 
matrix. 


4. Variable Relationships 

In the preceding text there have been two sets of relationships established. In 
Selection of Eigenvalues and Eigenvectors (Chapter III.A.2) a relationship between the 
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elements a, P and y of the scaled covariance matrix Cand the eigenvalues 5^. , s_ was 
established. A second relationship between the eigenvalues and the data parameters a, b, 0 
and was detailed in Error Ellipse Relationship to the Bivariate Normal Distribution 

(Chapter II.A). The following section uses these relationships to develop C in terms of a, 
b, 0 and x^- 

The eigenvector from the positive root contains information relative to the 
direction of the semi major axis of the ellipse. Of the two positive root solutions, it was 
determined that the second was the more numerically stable. Therefore, the x, y 

components of the eigenvector can be used to determine 0 as illustrated in Figure 3.3. 



From Equation 3.64 comes the following 

X(2) y 

xf? e,-p 

The semi-axes for the ellipse corresponding to a probability level p are related to the 
eigenvalues and s _. By Equations 2.22 and 2.23 we see that 
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(3.105) 

and 


1 - 

b = yjx^S. => - ^2 

(3.106) 

where x‘ « A® two-degree of freedom value of corresponding to probability 
Manipulating Equation 3.104 we find 

level p. 

-P = ycot0 

(3.107) 

which leads to the formula 


P=s^ -ycoXO. 

(3.108) 

From the Identity Equation 3.41 it follows that 


a = s^ +s_-p 

(3.109) 

= e^+s_-s^+ycoXd 

(3.110) 

= s_ +y cot 9 . 

(3.111) 

Substituting Equation 3.108 and Equation 3.111 into the Identity Equation 3.42 we find the 
second order polynomial 

s_ = ( 5 , + X cot 0)^s+ ~ Y cot 9 ) y 

(3.112) 

= s^e_ -ye_ GoXd+ye^ coXd-y^ cot^ d-y^ 

(3.113) 


This can be rewritten as 
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- (cot^ ^ +1) + cot $-e. cot ^) = 0 

(3.114) 

or 


/^{cot^ 0+l) + y{s. cot9-s^cot6) = 0. 

(3.115) 

Y {;'[cot^ ^ +1] + [ 5 . cot 9 -cot 0 ]} = 0 

(3.116) 

The non-zero solution of this equation is 


^^.cot^-£_ cot^ 
cot^^-hl 

(3.117) 

cot ^ 5 + 

1 

(3.118) 

sin^ 9 


= sin^ ^cot^fi; ~ ^-) • 

(3.119) 

Seen in a more convenient form throu^ application of the trigonometry 

identities as 

Y - sin^cos^fi; - ei). 

(3.120) 

Recalling Equations 3.105 looking for. Back and 3.106 we see that Equation 3.120 is y in 

terms of a, b, 9 and , exactly what we were substituting for a and 3 we 

; have 

a = ^ cot ^ 

(3.121) 

= g_ -i-sin^cos^fi'+ -5_)cot0 

(3.122) 

= £■_ -f- cos^ 9s^ - cos^ 9e_ 

(3.123) 


leading to the useful form 
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a = £■_ + cos^ d{s^ " • 


(3.124) 


For y0we get 

-ycotO (3.125) 

= s^- sin - e_)cotd (3.126) 

leading to the form we want 

P= s^- cos^ 6{s^ “ ^-) • (3.127) 

Equations 3.120, 3.124 and 3.127 establish a direct relationship between the 
covariance matrix and the data parameters a, b, 0 and . 


B. N+1 COMBINATORIAL EFFECTS 

The ability to combine multiple error ellipses, yielding the position vector and the 
associated covariance matrix is provided througji Equation 2.51. This section is concerned 
with the effect of including an additional error ellipse to the existing solution, i.e. the N+1 
case. In order to more clearly describe this problem, let us look at an example containing 
three ellipses. The question posed is will the solution for the case when all three 
observations are taken at once, as depicted in Figure 3.4, 
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be equal to the solution for the case in Figure 3.5, where two ellipsoids are initially 
processed, then a third is combined to that solution? 



Figure 3.5 COMBINATION EXAMPLE (b) 


Formally stated in Equation 3.128, we need to prove that die probability of the 
solution vector x given N+1 observations is equivalent to the probability of x given N 
observations given an additional observation. 

p(p(xloi...o^)|oj^+i)= P(xloi...o^^i) (3.128) 

The optimal estimate for the combination of N observed error eUipses was derived in 
Chapter II. Substituting the scaled covariance matrix C into Equation 2.51 we see that 
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f N \-V M _ '\ _ ^ ' 

Z«r’ Z«r'o, =«?; sc-'., 

/ V ,=1 / V ,=1 > 


W=1 


where 


c;= Z«r‘o, ^ er=Z«;' 

V i=i / /=i 


(3.129) 


(3.130) 


€7' is a new variable btroduced to differentiate between the previously summed and the yet 
to be summed N+1 cases. The scaled covariance matrix for the solution of N+1 
observations is 


=(cr +«».,) 


(3.131) 


^N+\ =(^tf ' +®v+i) 


(3.132) 


Now lef s work throu^ the N-h 1 solution 


Xv+i=(z®r‘] IS®r*o,- 

v,=i y vi=i ' 


(3.133) 


(S«r'o,+«i'.,o„.,j ( 3 . 134 ) 


= (e;-’ (3.135) 

V 1=1 


) N ^ ^ 

*V+1 “ +®v+l®v+l 

i=l 


(3.136) 


{c^ ’ +®v+lhv+l ~ ■'■®V+1®V+1 

' 1=1 


(3.137) 
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where = I 




+ ®//+l®JV+l 


(3.138) 


*Ar+l ~ (^AT ■•■®Af+l) I^AT *iV ®Af+l®Af+l] (3.139) 




•N+\\ 


I«r'o, 


' 1=1 


(3.140) 


Which agrees with Equation 3.128 proving the N+1 case is valid. This proof is critical 
from a software engineering point of view. Consider the N+1 case not valid. A piece of 
software wanting to keep a track updated over some amount of time, would need to store 
the information from all prior observations in order to update it’s location and associated 
error ellipse. Since the N+1 case is valid, only the most recent solution needs to be 
retained. Allowing all previous fixes to be discarded, freeing memory for other tasking. 


C. ELUPSE ROTATION 

A change of axes will be required in order to place the error ellipses that are to be 
combined in an approximately orthogonal coordinate system. In the latitude/longitude 
coordinate system, as latitude increases the coordinate system becomes progressively 
distorted. The closest approximation to an orthogonal coordinate system within the 
latitude/longitude coordinate system is found by transforming the ellipse center to a 
coordinate system where <j>=0, X=0. Where (j> is latitude and X longitude. This coordinate 
system will be referred to as the prime or coordinate system. 

The process will begin by defining the coordinate system axes in terms of 

(|> and A, with respect to the original ex,ey,e 2 coordinate system. Figure 3.6 depicts all 
axes and angles relevant to the transformation between the two coordinate systems. 
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Figure 3.6 ROTATION COORDINATE SYSTEM 


With the aid of Figure 3.6 the axis of the prime coordinate system and 

are defined as follows 


= cos(X)cos((l>)ex + sin(X)cos(<l>)ey + sin(<j))e^ (3.141) 

= cos(7i + X)cos^^- + sin(X + 7t)cos^Y“ (3.142) 

= -cos(/l)sin(^)ejj -sin(A)sin(<^)ey +cos(^^)e^ 


(3.143) 


As a simple check we verify that and are orthogonal by confirming their dot product 
is zero. 


.:_2 


- 61 ^ = - cos^ (X) sin((|») cos((j>) - sin^ (X) sin(<|)) cos(<l>) + sin(<j>) cos((j)) (3. 1 44) 


-[cos^ {X) + sin^ (X)] sin(^^) co + sin{^) cos(^) 


( 3 . 145 ) 





The third axis is easily found by the cross-product of the two previous directional vectors 




(3.146) 


= -sin(X)cos(>.)sin((l))cos(<l>)e^ +cos(>.)sin^((|))e^ +sin(X)cos(X)sin((l))cos(<|))4 


- sin(X,) sin^ ((|))ejc + cos(X) cos^ (<|))e^ - sin(X) cos^ (<t))e;c (3.147) 

= - sin(A,)ej. + cos(X)e^ (3.148) 

To ensure thorou^ness, the remaining combinations of axes are checked to be orthogonal 
e^-e^=- sin(X)cos(X,)cos(<l») + sin(X)cos(X,)cos((t>) = 0 (3.149) 

eQ-e^= sin(X,) cos(X,) sin(^) - sin(X,) cos(X) sin(<l)) = 0 (3.150) 


The unit vectors in the coordinate system can be express as a relationship 

of the directional cosines relative to the e^,e-f^,ei^ coordinate system. In unit terms 

4 ■ 4' = cos(4,4')»which leads to the following 


■ hh + (4 ■ ■ kh 

(3.151) 

cos(X.)cos(<|>)e^ -sin(?i)eT, -cos(X,)sin((|))e^ 

(3.152) 

(khh ■"’'V'n 

(3.153) 

sin(X)cos(<t>)e^ -i-cos(X)e^ - sin(X)sin(<l))e^ 

(3.154) 

= {e^ • e^)e^ + 

(3.155) 
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( 3 . 156 ) 


= sin(<|))e^ + cos((l))e(^ 




If we have coordinates x = 


y 

\zj 


in the coordinate system, the coordinates of the 

same point in the coordinate system are 




/*. A. 

< 

< 

r 

< 

< 


V 


e • e 

e •€ 

V y 

Cr •e 

C y 


T 

\zj 


a 

o 

o 

o 

sin(X,) cos(<t>) 

sin((t))^ 

rx\ 

- sin(X,) 

O 

o 

0 

y 

cos(X.) sin(<|)) 

- sin(X.) sin(<t)) 

cos(<t))y 

\z) 


The inverse transform is 


\ 


^AA A A 


X 


> 

>> 

> 

y 


^y'^r, ^y-^C 

rf 

\z) 


Jrh 

Id 


^cos(X,) cos((}>) - sin(A,) 
sin(A,)cos((|)) cos(X) 

sin((j)) 0 


cos(X,)sin(<|>)' 

V 

sin(X,)sin((j)) 

11 

a 

o 

Id 


(3.157) 


(3.158) 


(3.159) 


(3.160) 


If we let ({>i , Xi denote the latitude and longitude of a point in the coordinate 

system. Then x, can easily be expressed in terms (j>i and 


X,- = 

yi 


cos(</^.) 

C0s(/l5P 

lsin(A,) 


U,J 


^ sin 

M J 


(3.161) 
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The coordinates of the point in the ,^ti coordinate system are 




Ti 

f 


V 


cos(2)cos(^) 
- sin(A) 

- cos(/l) sin(^) 


sin(A)cos(^) 

cos(A) 

- sin(A) sin(^) 


sin(^)^ 

0 

cos(^)j 


^CO CO 

cos(^)sin(>^j) 

sin (< Z >,) 


(3.162) 


If denote the latitude and longitude of this point in the coordinate 

system, then 





^COs(^)cOs(>^5)^ 



= 

cos(<ij)sin(/l,) 




^ sin(^) J 


(3.163) 


Equation 3.163 can be used to derive the latitude and longitude of the point in the prime 
coordinate system. Table 3.1 shows the range or value of the longitude of the point after it 
is rotated. 


=sin ^(C/) 

2 2 

(3.164) ' 

I 

(1 

-7t < Xj <7t 

(3.165) 


Table 3.1 Rgcages of Prime Longitude 


If =0,11,. =0 

If < 0,T1, >0 —<Xi<Tl 

If =0,11, >0 Xi=^ 

If <0,Tij <0 -n<Xi <-— 

If > 0,r|,- >0 0<Xj < — 

If <0,T1,- =0 X./ = 7t 

If >0,T|/ <0 -—<Xj <0 

2 

If >0,Tly =0 X, =0 

If =0,T1,- <0 
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In order to transpose the heading of an error ellipse between the two coordinate 
systems we must know the relationship between the semi-major direction vector and the 
north pole. The coordinates of the north pole in the ex,ey,e 2 coordinate system are 

(3-166) 



The coordinates of this point in the coordinate system are found when Equation 

3.158 is applied to the north pole from the ex,ey,e 2 coordinate system 





'' cos(A)cos(^) sin(A) cos(<^) sin(^)^ 



Vnp 

\^np) 


- sin(A) cos(A) 0 

cos(A) sin(^) - sin(A) sin(^) cos(^)> 

0 

Kh 


(3.167) 


^sin(<|))^ 

0 

Uos((l)); 


(3.168) 


Using Equation 3.168 we find the latitude and longitude of the north pole in the prime 
coordinate system 


^np = sin ^(cos(<|))) 

0<d < — 

> — rnp — 2 

(3.169) 

T . -if 0 1 
^ Uinltl)); 

f0if<t>>0 
[It if (|)< 0 

(3.170) 

D. ELLIPSE HEADING RELATIVE TO 

NORTH POLE 



The individual error ellipse headings relative to true north are changed during the 
axes transformation. Consider an ellipse at some arbitrary latitude, with a heading of true 
north. When the ellipse is transformed to the prime coordinate system, the heading is no 
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longer valid as depicted in Figure 3.7. This section will present the process which calculates 
this change in heading. 



The first step in this process is to determining each error ellipse heading relative to 
true north . Let two points have latitude and longitude and We want to find 

the heading from A to B at point A measured relative to true north. The following pseudo 
code addresses the four trivial cases; one of the points is located on either of the 
poles, TV . First a point located on the north pole 


^ <f>A = 2 ’ = ^ 

^ then = 0 

For the south pole 

-f then = 0 
If = -f > 4 ^ -f then ^ 

Now the case 


If A, =4 

If ^.4 ^ 

a>AB=^ 
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Else If 

and the final trivial case, = X^+Tt 

IfA, =4±;r 

If 

Else If <-<l>B 

If no trivial cases are valid apply Napier’s analogies. 

With the trivial cases out of the way we move on to examine the cases where 
Napier’s analogies will need to be applied. Refer to Figure 3.8 for notation definitions used 
in the following problems. First consider the case where the longitude of point B is greater 

than that of point A 


AX, = — X,^ 

AX, > 0 ; AX, < 7t 

jt , , tt 


N 


Va I 

1 


A 

^ 1 1 

S 

B 


Figure 3.8 HEADING BETWEEN POINTS B AND A 


The points involved in this problem lie upon the surface of the earth. Taking this into 
account we must look to spherical trigonometry for a solution. Napier’s analogies (Selby, 
1967) provide the relationships from which this problem can be solved. Employed below 

we find that 
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Equation 3.173 becomes 
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(3.177) 


(3.178) 


Therefore we can find and cOg^ throng the following two formulas 


(0^5 = tan 


-1 


, N «’{jl'I'.l 

cotl-AXj ^ 


'1)1 


sinl -\^A + 


S + tan'M 


\ sin| 

-♦«]) 

J~ 

cos 



(3.179) 


a BA =tan *-1 


cotl — AX 
2 




Sin 




i^-tan'N 


1 


sixiU^A - 

cotl ^ AX - 

cos(^-[<t>^ + 


(3.180) 


Putting these angles in terms of heading we find that the heading from point A to point B is 


Heading A-> B = ( 0 ab 

and the heading from B to A is the mirror image of A to B 

Heading BA = 2n-(S >ba 


(3.181) 


(3.182) 
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Figure 3.9 illustrates the case when the longitude of A is larger than the longitude of B 


AX = — Xq 

AX > 0; AX < TC 

< 1 >^ 


Figure 3.9 HEADING BETWEEN POINTS A AND B 

Here the same formulas for ©ab and ©ba apply. The heading from point A to B is 

Heading = (3.183) 

and the headbg from B to A is the mirror image of A to B 

Heading 5 ^4 = © 34 (3.184) 

Now that we understand the how to obtain the heading between two points the 
next step is to apply this theory to our two coordinate systems. Let ©,• denote the heading 

from the center of each ellipse (j);, X,- to the north pole ^ in the ^Tji^-coordinate system. 

To find ©/, let $/ ,Xj = point A and = point B and execute the pseudo code in the 

beginning of this section. If 0; is the heading of the ellipse at the point (})/ ,Xj in the xyz- 

coordinate system, then the heading of the ellipse at the point <|)j, X,- in the ^rj^-coordinate 
system is the sum of the original heading and the offset caused by rotation 

e,=®/+0/ (3.185) 
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IV. COMBINATION ALGORITHM 


A. PROCESS OUTLINE 

The Combination Algorithm will apply the theory and concepts from Chapters II 
and III to the real world environment. The input parameters for the combination 
algorithm will consist of a set of error ellipses with center , Aj (/=!, ... , N), semi-major 
axis A, semi-minor axis B, , heading Bt and probability P; (=> ) • Th® algorithm outlined 

below will process the N input ellipses and produce a solution ellipse that represents an area 
that has the specified probability that the target lie within it’s bounds. Following is the 
Combination algorithm presented in a step by step manner. 


1. Input a set of error ellipses with center , A, (/=1, ..., N), semi-major axis A , 
semi-minor axis B, , heading 9i and probability p, (=> xf ) » where , A/ refer to the latitude 

and longitude respectively. The semi-axes A/ and B/ are in an appropriate unit of length, Bt 
in degrees and the desired probability level p refers to the confidence that the target lie 
within the generated solution ellipse. 

2. Project semi-axes A and B, onto a plane to obtain a, and b/, where 




a, = r, tan] 


KrJ 


(4.1) 


b,= 



(4.2) 


Rationalization and development of Equation 4.1 and Equation 4.2 is presented in 
Chapter III.A.1, Scaled Covariance Matrix. 

3. Find the eigenvalues corresponding to the semi-major axes a,, semi-minor axes b, 
and the probability level p of each ellipse, where 
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(4.3) 


^+(0 ■" 




(4.4) 


Equation 4.3 and Equation 4.4 find their bases in Chapter II.A, Relationship of the 
Error Ellipse to a Bivariate Normal Distribution. The values of p for the distributed 
function are given in Table 2.1. 

4. Find the first approximation for the combined location of the center of the 

ellipses 


1 ^ 

^ i=\ 

(4.5) 

1 ^ 

= ■^Scos(^)sin(A) 

ly i=\ 

(4.6) 



.(4-7) 


I' 


«!> = 

. -1 z 
sm — 

(4.8) 




X 



(4.9) 


The first part of this step finds the simple average of the (x,y, 2 ) components of the 
individual error ellipses. Note that the information within the individual observation 


vectors x,, where x,- = 




T, 

\z,) 


, should be retained for future use. These x, vectors will be 


needed for Step 5. The second part of this step converts the averaged (X/^<,^() components, 
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producing the latitude (j> and longitude X which is an approximation of die center of all the 
observed error ellipses. 

5. Transform the location of the center of each error ellipse from xyz-coordinates 
to (^74~coordinates (prime coordinates). 


^ cos(A) cos(^) 

sin(A) cos(^) 

sin(^)^ 

cos(^)cos(a,P 


- sin(A) 

cos(/l) 


0 

r, cos(«^)sin(;ii) 

(4.10) 

cos(A) sin(^) 

- sin(A) sin(^) 

a 

o 


reSin(<!^) V 


^ cos(A) cos(^) 

O 

o 

sin(^)^ 




- sin(A) 

cos(A) 

0 

yt 


(4.11) 

V- cos(2) sin(^) 

- sin(A) sin(^^] 

O 

o 

U,J 




<t>i = sin ’ 

VrJ 




(4.12) 


X^ = tan"’ 





(4.13) 


Note that x, calculated in Step 4 does not need to be recalculated to obtain 4, , 


where 


" 4 ^ 


The transformation matrix used in Equation 4.10 is based on Equation 


3.158 developed in Chapter III.C, Ellipse Rotation. During this step the first approximation 
coordinates ^,/l (Step 4) are the focus upon which the coordinate system is rotated as 
illustrated in Figure 4.1. The end product being the spatial relationship of each observation 
^,X to the first approximation ^is maintained, but transformed to a coordinate system 
where ^ = 0,A = 0. This transformation is required compensate for the difference in arc- 
length due to a fixed change in longitude as the latitude increases from zero. 
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First Approximation 
Coordinates 



6. Find the coordinates of the north pole in ^//^-coordinates 



(4.14) 



0 if <t>>Q 
n if ^ < 0 


(4.15) 


Location of the north pole in prime coordinates is needed to solve for die new 
ellipse heading after the coordinate rotation. Figure 4.2 shows the geometry involved in the 
pole rotation. Equations 4.14 and 4.15 were obtained form Equations 3.169 and 3.170. 


54 





Figure 4.2 NORTH POLE ROTATION 

7. Find the heading from the center of each error ellipse ^to the north 
pole ip in ^T/C-coordinates. 

The process for this step is detailed in Chapter III.D, Ellipse Heading Relative to 
North Pole. In execution, the first move is to examine the trivial cases A, =2^^ , 

^ = 2„p + 7 r, and the case where one of the points is located on the north pole. In the 
likely event that none of the trivial case match, Napier’s analogies must be utilized by 
employing Equation 3.179. 

8..Find the heading of the error ellipses in It/^^coordinates 

^ = (4.16) 

The heading of the individual error ellipses in the prime coordinate system is the 
sum of the heading of each error ellipse S, (Step 7) and ellipse’s heading in xyz-coordinates 

(Step 1). 

9. Find the scaled covariance matrix corresponding to each error ellipse 
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II 


(4.17) 



(4.18) 

a. =£_(,)+cos'^ 

(4.19) 

1 

o 

o 

1 

il 

(4.20) 

=sin^ 

1 

(4.21) 


Along with the scaled covariance matrix C,. , Step 10 will require the individual 

ellipse scaled covariance matrix inverse C:’ also be calculated. The elements of these 
matrices are derived in Chapter III.A.3, Variable Relationships. 

10. Find the 2x2 scaled covariance matrix for the sum of the error ellipses 


c-=sc-' = 


1=1 


'a 

<v r) 


( 4 . 22 ) 



ar-u \-o 



a 

r 


(4.23) 


Equation 4.22 is a result of the methodology outlined in Chapter II.C. The 2x2 

scaled covariance matrix for the sum of the error elUpses C produced in this step 
represents the solution ellipse that has the specified probabiUty that the target lie within it’s 
bounds. The location and heading of this solution ellipse has yet to be determined. 

11. Convert the 2x7 scaled covariance matrix for the sum of the ellipses to the 3x3 
matrix 
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1 


S"' 


r 


0 

1 

r^cos^j 


(4.24) 






(4.25) 



r 



\ 0 


^ix 0 

0 o; 




di^ 


dX 

dh 

dr] 

drj 

drj 


dx 

~dh 




\d^ 

dX 

dh) 




(4.26) 


(4.27) 


(4.28) 


There are four tasks to accomplish prior to producing the 3x3 scaled covariance 
matrix for the sum of the ellipses ^ . First, the inverse of the scaling matrix defined in 

Equation 3.7 must be generated. Second, using the 2x2 scaled covariance matrix C firom 
Step 10, calculate the unsealed covariance matrix in terms of latitude and longitude. Next, 
populate the 3x3 unsealed covariance matrix. This matrix C^~ ^ is in terms of latitude, 

longitude and hei^t. Lastly, the 3x3 scaling matrix T needs to be produced. T contains 
the partial derivatives of ^,7 and ^ with respect to A and h (Equations 3.19 thru 3.27). 
With those tasks complete can be calculated. This conversion process is described 

in detail in Chapter III.A.2.a 2x2 to 3x3 Covariance Matrix Conversion. 
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12.Find the wei^ted average location for the center of the error ellipse in ^rj^— 
coordinates 


i^] 

N 





Id 

J=1 



^ = sin 


-1 


r 




X 


= tan ‘ 



(4.29) 

(4.30) 

(4.31) 


The vector 4 represents the location of the weig^ited average of the N 
observations. The theory behind Equation 4.29 is presented in Chapter II.C, Methodology. 
The new latitude (ft and longitude A are in essence refinements of the rotated first 
approximation. The level of refinement will be monitored once the solution is transformed 
back to the xyz-coordinate system. If convergence upon the true solution is not obtained, 
an iterative process will begin. 

13. Find the eigenvalues of the combined covariance matrix. 

+1 [(a - + 4;^" (4.32) 




(4.33) 


where a,P and y are obtained from Equation 4.23. 

The eigenvalues produced in this step are the key to extracting information about 
the semi-major axis, semi-minor axis and heading of the covariance matrix. Equation 4.32 
and Equation 4.33 supply the real solutions to the characteristic polynomial of the 
covariance matrix discussed in Chapter III.A.3.a, Eigenvalues. 
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14. Find the heading of the error ellipse in the ^ 77 ^-coordbate system 


e 


= tan 



Vs^-P) 


(4.34) 


The error ellipse heading is derived from the ^ and rj components of the 
eigenvectors found in Section III.A.3.C, Eigenvectors. The development of Equation 4.34 
can be found in Chapter III.A.4, Variable Relationships. 

15. Find the heading a from (ff, X to (f>„p, X„p. 


The same process found in Step 7 is utilized here. 

16. Find the heading of the error ellipse in the xyz-coordinate system 

B = 9 - CO (4.35) 


The heading of the error ellipse in the xyz-coordinate system is found subtracting 
the h ead ing offset due to rotation from the ellipses heading in the prime coordinate system. 
Note that Equation 4.35 agrees with Equation 4.16 from Step 8. 

17. Find the coordinates of the weighted average location for the center of the 
error ellipses in xyz-coordinates 


X 


0 

0 

0 

0 

- sin(A) 

- cos(A)sin(^)^ 


y 

=: 

sin(A)cos(^) 

cos(A) 

- sin(/l)sin(^) 

V 

UJ 


1 sin(<i>) 

0 

cos(^) ; 

Id 


(3.36) 



(4.37) 


A = tan ‘I t 


(4.38) 
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The transformation from prime coordinates to the xyz-coordinate system is 
developed in Chapter III.C, Ellipse Rotation. This step is the inverse of Step 5. 

18. Check for convergence 


If \(I)-<I,\<K^ and 

else set 


<Kx then 


goto Step 19 
x = x ; y=y ; z=z 


(f) ^ (j> X — X 

and goto Step 5 


In some cases the difference between the first approximation and the 
weighted average location my be significant. By setting values for and that are small 

an iterative process is established to fine tune the solution. 

19. Calculate the semi-major and semi-minor axes of the combined error ellipse for 
the derived probability level 


a = 

(4.39) 


(4.40) 


The semi-major and semi-minor axes of the combined error ellipse are found by 
utilizing the eigenvalues calculated in Step 14 and a selected chi-square number for Table 
II.l. Chapter II.A, Relationship of the Error Ellipse to A Bivariate Normal Distribution 
gives the origins of Equation 4.39 and Equation 4.40. 

20. Project the semi-major and semi-minor axes from the plane to a spherical earth 


A = 


r, tan 


-1 

UJ 


(4.41) 
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(4.42) 


B = r. tan"' 


-1 

UeJ 


The last step in the Combination Algorithm places solution ellipse back onto the 
spherical earth. This reverses Step 2 and supplies the last parameters of the solution ellipse 

which consists of latitude (f>, longitude X, semi-major axes A, semi-minor axes B, and 
heading 9. 


B. PROCESS FLOW CHART 

Figure 4.3 provides a condensed view of the algorithm. For details of specific steps 
refer to Section A of this chapter. 



Figure 4.3 PROCESS FLOW CHART 
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V. SUMMARY & CONCLUSIONS 


A. AREAS FOR FUTURE STUDY 

One of the key assumptions made in the development of this algorithm was that the 
error ellipse data is derived from a bivariate normal distribution. Several transformations 
were made throu^out the derivation to transform spherical data to a form that could be 
analyzed using a bivariate normal distribution. When the error ellipse is small (a << r_e) the 
probability distribution is highly localized and these transformations lead to a good 
approximation of a spherical probability distribution. The algorithm developed here will 
therefore be applicable to most error ellipses obtained from operational systems. Strictly 
speaking, however, we should have used a probability distribution more suited to spherical 
data, the Kent distribution (Fisher, Lewis, Embletonl987). While this distribution is more 
difficult to work with, it would be valid for error ellipses of all sizes. This issue should be 
examined more fully in the future. 

There are three remaining areas relevant to this thesis that are available for future 
research. The first and foremost is to design and implement a test and validation program. 
The validation program would provide a means to quantitatively verify the theory put forth 
in this research. The project would entail producing computer code to implement the 
Combination Algorithm. Real world data sets are available to compare against output that 
would be generated by the program. 

Broadening the Combination Algorithm to include multiple source data is the next 
potential research area. The ability to correlate observations from multiple platforms would 
greatly enhance our national warfi^ting capability. Decision makers would no longer need 
to mentally extrapolate a solution from numerous observations of die same target. The 
generation of this composite geolocation would reduce total decision cycle time, ultimately 
leading to increased force effectiveness. With an effective combination algorithm in place 
the need to transmit multiple geolocations of the same target would be eliminated, reducing 
fleet bandwidth requirements. By combining observations of the various platforms both 
bandwidth and data storage requirements are reduced. The reduction would be a result of 
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two factors: (1) collection platforms will combine same source observations reducing the 
amount of data to be transmitted, (2) once the collection platform solutions are received 
and combined into a composite solution the data received from the collection platforms can 
be discarded. The problem of including multiple source data could be approached by 
expanding the observation vector to include system dependent biases such as atmospheric 
delays and refraction. 

The third area to consider for future research is the correlation of observations 
obtained of a moving target. Expanding the domain of possible targets from fixed site to 
include mobile platforms would significantly reduce the information presented to a decision 
maker. By expanding the number of platforms available for data reduction throu^ 
combination, bandwidth and data storage resources are conserved. The inclusion of 
moving targets mi^t be accomplished by sampling and processing a discrete number of hits 
over a large number of observations. Those interested in pursuing any one of these area 
can contact NRaD Code 841 for information on current research. 

Areas of interest observed during initial attempts to implement the Combination 
Algorithm revolved around data structure design. Because the vast majority of operations 
to be performed in the algorithm involve matrix arithmetic, a language designed to for that 
purpose such as MATLAB is recommended. If C++ is to be used, it is recommended that 
an established toolbox of data objects and algorithms such as M++ by Dayd Software be 
utilized. Because of the large number of multi-dimensioned arrays required to code this 
algorithm a large memory model will be required. In the event that a smaller memory 
model must be used, the number of observations the program will be able to process will 
be severely restricted. To further reduce the strain on memory resources dynamic allocation 
of objects should be incorporated. 
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B. CONCLUSION 

This thesis has set the theoretical basis and developed an algorithm for optimally 
combining multiple small spherical error ellipses into a single error ellipse. This algorithm is 
applicable to most of the error ellipses produced by operational systems. When diis 
algorithm is validated it could have a major impact on how geolocation data is presented to 
the warfighter. In the future, as the electronic density of the battle field increases, the need 
for fusion of geolocation data will become increasingly apparent. This work can be 
considered the first step toward providing the warfi^ter with a fused multi-sensor 
geolocation solution. 
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